The pseudospectral method : Accurate representation in elastic wave calculations

When finite-difference methods are used to solve the elastic wave equation in a discontinuous medium, the error has two dominant components. Dispersive errors lead to artificial wave trains. Errors from interfaces lead to circular wavefronts emanating from each location where the interface appears *‘jagged” to the rectangular grid. The pseudospectral method can be viewed as the limit of finite d&lerences with infinite order of accuracy. With this method, dispersive errors are essentially eliminated. The mappings introduced in this paper also eliminate the other dominant error source. Test calculations contirm that these mappings significantly enhance the already highly competitive pseudospectral method with only a very small additional cost. Although the mapping method is described here in connection with the pseudospectral method. it can also be used with high-order finite-difference approximations. INTRODUCTION The pseudospectral method was introduced in the early 1970s (Kreiss and Oliger. 1972; Orszag, 1972; Fornberg, 1975). It has since gained wide acceptance in a variety of areas, including weather forecasting, nonlinear waves, and turbulence modeling The pseudospectral method did not impa~~m the field of seismic modeling until about 10 years after its introduction into other areas (Kosloff and Haysal, 1982; Johnson, 19X4), most likely because, initially, neither nonperiodic domains nor discontinuous interfaces appeared practical with the method. However. these difficulties are gradually being overcome. Fornberg (1987) discusses the basic features of the pseudospectral method and compares the method to finite-difference methods for the 2-D elastic wave equation. In the case of of interfaces smoothly varying coefficients and with interfaces aligned with the grid, it was found that the required grid spacings in each space dimension satisfied ratios of 16:4 : 1 for the pseudospectral method, fourth-order finite differences, and secondorder finite differences. These ratios were derived theoretically from an analysis of dispersion errors in homogeneous media and were confirmed in more general test calculations. Issues relevant to practical implementations and not directly addressed in my earlier work include the following: (i) For the extremely coarse grids which are sufficient to carry traveling pulse solutions in the pseudospectral method. some technique must be found to specify interfact locations more accurately than to the nearest grid point. The roughness in the discretization of a smooth interface not aligned with the grid acts as a source of unphysical noise. (ii) The performance of the pseudospectral method for surface and interfacial waves must be examined. (iii) Nonreflecting boundary conditions which are compatible with the periodicity in space implicit in the pseudospectral method must be devised. (iv) time integration techniques must be implemented in a way which minimizes computer time and storage. Adequate, but most likely not optimal. solutions to the last two problems have been presented elsewhere. One way to address problem (iii) is to extend the computational region with absorbing boundaries. Cerjan et al. (1985), KoslofS and Kosioil‘ ( 1986). and Loewenthal et al. ( 1987) present test results of such implementations. Regarding problem (iv), note that straightforward central ditrerencing in time (leap-frog) is only accurate to second order and requires two time levels of storage. assuming the 2-D equations are formulated as five coupled first-order equations. Higher order accuracies can be achieved by the “method of lines” approach. Splitting methods (e.g., Bayliss et al., 1986) are limited to second-order accuracy, but can reduce the storage requirement to one time level. Manuscript received by the Edltor March 9, 1987; revised manuscript received September 29, 1987 *Exxon Research and Enrmeering Comoanv. Annandale NJ OXROI. (‘ 11)8R Society of Explor&n Ge%phys&t~. All rights reserved. 625 Pseudospectral Interface Representation 627 J’(<, rl + 1) = Yk 11) + 1. Table 1. Number of arithmetic operations for each time step with the different methods.